##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--gene_file"), type = "character") ,
    make_option(c("--smg_file"), type = "character") ,
    make_option(c("--cgc_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.tsv",sep="")
    info_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv",sep="")
    smg_file <- paste(work_dir,"/public_ref/importTantGene.list",sep="")
    cgc_file <- paste(work_dir,"/public_ref/cancer_gene_census.list",sep="")
    gene_file <- "~/20220915_gastric_multiple/dna_combinePublic/Pyclone/selectGCClone/GCClone_gene.reord.tsv"
    images_path <- "~/20220915_gastric_multiple/dna_combinePublic/Pyclone/selectGCClone"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

info_file <- opt$info_file
mut_rate_gene_file <- opt$mut_rate_gene_file
images_path <- opt$images_path
gene_file <- opt$gene_file
cgc_file <- opt$cgc_file
smg_file <- opt$smg_file

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")
col_im_gc <- col[c(1,4)]
col_igc_dgc <- col[c(2,3)]

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_smg <- data.frame(fread(gene_file))
dat_info <- data.frame(fread(info_file))
smg <- dat_smg$gene

dat_report <- data.frame(fread(smg_file))
dat_cgc <- data.frame(fread(cgc_file))

###########################################################################################

tsg <- dat_cgc[grep( "TSG" , dat_cgc$Role.in.Cancer),1]
oncog <- dat_cgc[grep( "oncogene" , dat_cgc$Role.in.Cancer),1]
tsg_oncog <- tsg[tsg %in% oncog]

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol %in% smg)

dat_plot <- rbind( dat_mutRateGene )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IM" , "GC" , "IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 

###########################################################################################

dat_plot <- subset( dat_plot , From == "All" )

im_num <- unique(subset(dat_plot , Class=="IM")$SampleNum)
gc_num <- unique(subset(dat_plot , Class=="GC")$SampleNum)
igc_num <- unique(subset(dat_plot , Class=="IGC")$SampleNum)
dgc_num <- unique(subset(dat_plot , Class=="DGC")$SampleNum)

###########################################################################################
## 注释样本的分子亚型
dat_info <- subset( dat_info , Class != "IM" )
dat_info <- unique(dat_info[,c("Normal" , "TCGA_Class")])
rownames(dat_info) <- dat_info$Normal

###########################################################################################
## 计算P值
result <- c()

for(geneN in unique(dat_plot$Hugo_Symbol)){

    print(geneN)

    tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("GC") )
    tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("IM") )

    tmp_igc <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("IGC") )
    tmp_dgc <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("DGC") )

    if(nrow(tmp_1)==0){
        tmp_1 <- tmp_2
        tmp_1$MutNum <- 0
        tmp_1$MutRate <- 0
        tmp_1$value_text <- ""
    }

    if(nrow(tmp_2)==0){
        tmp_2 <- tmp_1
        tmp_2$MutNum <- 0
        tmp_2$MutRate <- 0
        tmp_2$value_text <- ""
    }

    if(nrow(tmp_igc)==0){
        tmp_igc <- tmp_1
        tmp_igc$MutNum <- 0
        tmp_igc$MutRate <- 0
        tmp_igc$value_text <- ""
    }
    if(nrow(tmp_dgc)==0){
        tmp_dgc <- tmp_1
        tmp_dgc$MutNum <- 0
        tmp_dgc$MutRate <- 0
        tmp_dgc$value_text <- ""
    }

    tmp <- rbind( tmp_1 , tmp_2 )

    tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
    p <- fisher.test(tmp_fisher)$p.value

    tmp <- data.frame( gene = geneN , MutRate_IM = tmp_2$MutRate , MutRate_GC = tmp_1$MutRate , 
        MutRate_IGC = tmp_igc$MutRate , MutRate_DGC = tmp_dgc$MutRate , p_IM_GC = p
        )

    result <- rbind( result , tmp )
}

dat_out <- merge( result , dat_smg , by = "gene" , all.x = T , all.y = T )
dat_out$Filter3 <- ifelse( dat_out$p_IM_GC <= 0.05 & dat_out$MutRate_IM >  dat_out$MutRate_GC , "IM_HighMutRate" , dat_out$Filter2 )

###########################################################################################
## 注释样本的分子亚型
dat_out$Molecular.subtype <- ""
for(geneN in unique(dat_out$gene)){

    print(geneN)
    tmp <- subset( dat_out , gene == geneN)
    normals <- unlist(strsplit(tmp$Normal , ","))
    dat_out[which(dat_out$gene == geneN),"Molecular.subtype"] <- paste0(dat_info[normals,]$TCGA_Class , collapse = ",")
}

col_order <- c("gene" , 
    "MutRate_IM" , "MutRate_GC" , "MutRate_IGC" , "MutRate_DGC" , "p_IM_GC" ,
    "Class" ,
    "CloneMutNum" , "Normal" , "Molecular.subtype" , 
    "AddtionalCloneSampleNum" , "AddtionalCloneSample" ,
    "Filter1" , "Filter2" , "Filter3"
    )

dat_out <- dat_out[,col_order]

###########################################################################################
## 标记基因为SMG还是CGC
dat_out[dat_out$gene %in% tsg & dat_out$Filter3 == "SMG" & dat_out$Filter3 != "IM_HighMutRate"  , "Filter3" ] <- "TSG"
dat_out[dat_out$gene %in% oncog & dat_out$Filter3 == "SMG" & dat_out$Filter3 != "IM_HighMutRate"  , "Filter3" ] <- "Oncogene"
dat_out[dat_out$gene %in% tsg_oncog & dat_out$Filter3 == "SMG" & dat_out$Filter3 != "IM_HighMutRate"  , "Filter3" ] <- "TSG_Oncogene"
dat_out[dat_out$gene %in% dat_report$Gene_Symbol & dat_out$Filter3 != "IM_HighMutRate"  , "Filter3" ] <- "SMG"

dat_out$Filter_use <- ifelse( dat_out$Filter3 %in% c("Clone_Independent & NoSMG_1Clone & AddtionalCloneSample" , 
    "Clone_Independent & NoSMG_2Clone" , "Oncogene" , "SMG" , "TSG") , "TRUE" , "False" )

images_name <- paste0(images_path , "/GCClone_gene.reord.addMutRate.tsv")
write.table( dat_out , images_name , row.names = F , sep = "\t" , quote = F )
